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Abstract 

We develop a fast algorithm to construct the robustness degradation function, which describes 
quantitatively the relationship between the proportion of systems guaranteeing the robustness 
requirement and the radius of the uncertainty set. This function can be applied to predict whether 
a controller design based on an inexact mathematical model will perform satisfactorily when 
implemented on the true system. 

1 Introduction 

In recent years, there has been growing interest in the development of probabilistic methods for 
robustness analysis and design problems aimed at overcoming the computational complexity and the 
conservatism issue of the deterministic worst-case framework [T5J El El El El El El SI El El El 
El [2T| [2^1 [151 E]- In the deterministic worst-case framework, one is interested in knowing if the 
robustness requirement is guaranteed for every value of the uncertainty. However, it should be borne 
in mind that the uncertainty set may include worst cases which never happen in reality. Instead of 
seeking the worst-case guarantee, it is sometimes "acceptable" that the robustness requirement is 
satisfied for most of the cases. It has been demonstrated that the proportion of systems guaranteeing 
the robustness requirement can be close to 1 even if the radii of the uncertainty set are much 
larger than the worst-case deterministic robustness margin [21 [13| (H El HHJ- Therefore, it is of 
practical importance to construct a function which describes quantitatively the relationship between 
the proportion of systems guaranteeing the robustness requirement and the radius of the uncertainty 



'Received by the editors June 5, 2002; accepted for publication (in revised form) August 6, 2003; published elec- 
tronically January 22, 2004. This research was supported in part by grants from NASA (NCC5-573) and LEQSF 
(NASA/LEQSF(2001-04)-01). 



1 



set. Such a function can serve as a guide for control engineers in evaluating the robustness of a 
control system once a controller design is completed. Such a function, referred to as a robustness 
degradation function, has been proposed by a number of researchers [21 [8]. For example, Barmish, 
Lagoa, and Tempo [2] have constructed a curve of the robustness margin amplification versus risk in 
a probabilistic setting. In a similar spirit, Calafiore, Dabbene, and Tempo [U [9] have constructed a 
probability degradation function in the context of real and complex parametric uncertainty. 

In this paper, allowing the robustness analysis to be performed in a distribution-free manner, 
we introduce the concept of proportion and adopt the assumption from the classical robust control 
framework that uncertainty is deterministic and bounded. It follows naturally that the robustness 
of a system can be reasonably measured by the ratio of the volume (Lebesgue measure) of the 
set of uncertainty guaranteeing the robustness requirement to the overall set of uncertainty |19| . 
Evaluation of such a measure of robustness requires generating samples with uniform distribution 
over uncertainty sets such as a spectral normal ball or an L ball. The difficulty of generating such 
samples has been successfully resolved in [S (9] . 

The conventional method for constructing the robustness function is to perform, independently, 
a certain number of simulations for each value of the uncertainty radius and then plot the function. 
Although such a curve can be applied to evaluate the robustness of the control system, it may be 
computationally expensive. This is especially true when many cycles of controller synthesis and 
robustness analysis are needed in the development of a high performance control system. Motivated 
by this situation, we focus on the machinery that can make the construction of such a function 
efficient. We have developed a sample reuse algorithm that allows the simulations to be conducted 
in an iterative manner. The idea is to start simulation from the larger uncertainty set and save 
appropriate evaluations of the robust requirement for the use of later simulations on the smaller 
uncertainty set. In this way the total number of simulations can be reduced significantly as compared 
to the conventional method. 

In addition to deriving our sample reuse algorithm from the worst-case deterministic framework, 
we show that the technique is also applicable when considering the random nature of the uncertainty. 
In such cases, the worst-case properties of uniform distribution given in the pioneering work [5J [21 Q] 
allow our algorithm to be applied to efficiently solve a wide variety of robustness analysis problems. 
In particular, the radial truncation theory [2] can be applied to robustness analysis problems with 
uncertainty bounding sets defined as spectral norm balls and l p balls. 

The organization of the paper is as follows. Section 2 gives the problem formulation. Section 3 
presents our sample reuse algorithm. Section 4 is the performance analysis of the algorithm. Section 
5 applies the algorithm to examples. Section 6 shows the justification of the algorithm for the case 
of random uncertainties. Section 7 is the conclusion. The proofs of the theorems are included as an 
appendix. 



2 



2 Problem formulation 



We adopt the assumption, from the classical robust control framework, that the uncertainty is 
deterministic and bounded. We formulate a general robustness analysis problem as follows. 

Let P denote a robustness requirement. The definition of P can be a fairly complicated combi- 
nation of the following: 

• stability or Instability; 

• Hqo norm of the closed-loop transfer function; 

• time specifications such as overshoot, rise time, settling time, and steady state error. 

Let B(r) denote the set of uncertainties with size smaller than r. In applications, we are usually 
dealing with uncertainty sets such as the following: 

• l p ball B p {r) := {A € R n : ||A|| P < r} , where denotes the l p norm and p = 1, 2, . . . , oo. In 
particular, Boa{r) denotes a box. 

• Spectral norm ball B a (r) := {A £ A : <r(A) < r}, where cx(A) denotes the largest singular 
value of A. The class of allowable perturbations is 



where qi € F, i = l,...,s are scalar parameters with multiplicity n,...,r s and A,; E 

F niXm % i = l,...,c are possibly repeated full blocks. Here ¥ is either the complex field 
C or the real field R. 

• Homogeneous star-shaped bounding set Bu{r) := {r(A — Ao) + Ao : A G Q} , where Q C R n 
and Ao € Q (see [2] for a detailed illustration). 

Throughout this paper, B(r) refers to any type of uncertainty set described above. Define a function 
£(.) such that, for any X, 



A := {blockdiag[gi/, 



• • • , q s ir s , Ai 



...,A C ]} 



(1) 



£(X) := mm{r : X € B(r)} 



i.e., B(£(X)) includes X exactly in the boundary. By such definition, 




£{X) = a(X) 
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and 

£{X) = \\X\\ P 

in the context of a homogeneous star-shaped bounding set, spectral norm ball, and l p ball, respec- 
tively. 

To allow the robustness analysis to be performed in a distribution-free manner, we introduce the 
notion of proportion as follows. For any A € B(r) there is an associated system G(A). We define 
proportion as follows: 

vol({A 6 B(r) : The associated system G(A) guarantees P}) 

P r := utz( \\ 

vol(o(r)J 

with 

vol(S) := / dq, 
Jqes 

where the notion of dq is illustrated as follows: 



• (I): If q = [x rs ] nxrn is a real matrix in R nxm , then dq = Ylr=i ECli dx rs . 

• (II): If q = [x rs + jy rs ]nxm is a complex matrix in C nxm , then dq = Yl™=i W^=i{dx rs dy rs ). 

• (III): If q G A, i.e., q possesses a block structure defined by (P), then dq = (PIi=i ^Qi)(rii=i ^«)> 
where the notion of dg^ and dA j is defined by (I) and (II) . 

It follows that P(r) is a reasonable measure of the robustness of the system [51 [2D]. In the worst- 
case deterministic framework, we are interested only in knowing if P is guaranteed for every A. 
However, one should bear in mind that the uncertainty set in our model may include worst cases 
which never happen in reality. Thus, it would be "acceptable" in many applications if the robustness 
requirement P is satisfied for most of the cases. Hence, due to the inaccuracy of the model, we should 
also obtain the value of P(r) for uncertainty radius r which exceeds the deterministic robustness 
margin. 

Clearly, P(r) is deterministic in nature. However, we can resort to a probabilistic approach to 
evaluate P(r). To see this, one needs to observe that a random variable with uniform distribution 
over B(r), denoted by A", guarantees that 

Pr{A« e S} = V ° 1( ^ g y } 
vol(o(r)J 

for any S, and thus 

P(r) = Pr{The associated system G(A") guarantees P}. 



4 



It follows that a Monte Carlo method can be employed to estimate P(r) based on independently 

and identically distributed (i.i.d.) observations of A". 

It is interesting to know how the function P(r) degrades with respect to r when r increases from 
a to b, where b, a > 0. In a similar spirit, such a function has been proposed as a confidence 
degradation function in [2] and as a probability degradation function in [HI [9]. In this paper, we 
refer to the function P(.) as a robustness degradation function for the following reasons. First, we 
introduce the confidence interval for assessing the accuracy of the estimate of P(r). To be useful, 
every numerical method should be associated with an assessment for the accuracy of the estimate. 
Monte Carlo simulation is no exception. To avoid confusion, we reserve the notion of "confidence" 
for the purpose of interval estimation. Second, we introduce the concept of proportion for measuring 
robustness, which has no probabilistic content. Third, P(r) is a robustness measure and is usually 
decreasing with respect to r when P(r) is close to 1. 

To construct such a function of practical importance, the conventional way is to grid the interval 
[a, b] as a = p\ < p2 < ■ ■ ■ < pi = b and estimate P(/?i) by conducting N i.i.d. sampling experiments 
for each pi. In total, we need Nl samples. In the next section we show that the number of experiments 
can be significantly reduced. 

3 Sample reuse algorithm 

To improve efficiency, we shall make use of the following simple yet important observation. 

Let q* be an observation of a random variable with uniform distribution over B{r*) D <S(r) such 
that q* £ B(r). Then q* can also be viewed as an observation of a random variable with uniform 
distribution over B{r). 

In our algorithm, we flip the order of pi by defining 

n = pi+i-i 

for i = 1,2,... , I. Thus, the direction of simulation is backward. Our algorithm is described as 
follows. 

Sample Reuse Algorithm 

• Input: Sample size N, confidence parameter 5 6 (0, 1) and uncertainty radii r^, i = 1, 2, . . . , I. 

• Output: Proportion estimate Pj and the related confidence interval for i = 1, . . . ,1. In the 
following, run denotes the number of sampling experiments conducted at rj, and % denotes 
the number of observations guaranteeing P during the mn sampling experiments. 

• Step 1 (initialization). Let M = [mjj]; X 2 be a zero matrix. 
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• Step 2 (backward iteration). For i = 1 to i = I do the following: 

— Let r <— n. 

— While mil < -/V do the following: 

* Generate uniform sample q from £>(r). Evaluate the robustness requirement P for q. 

* Let m s i <— m s i + 1 for any s such that r > r s > £(q). 

* If robustness requirement P is satisfied for q, then let m S 2 <— m S 2 + 1 for any s such 
that r > r s > i(q). 

— Let Pj <— and construct the confidence interval of confidence level 100(1 — 5)%. 

It follows that q can be viewed as an observation of a random variable with uniform distribution 
over B(rj) if and only if r > r,- > £(q). Hence, if the robustness requirement P has been evaluated for 
£>(rj) at sample q, the result can be accepted without repeated evaluation of P for all B(rj) such that 
r > rj > i(q). Thus, sample reuse allows us to save both the sample generation and the evaluation 
o/P for the sample. It is also interesting to point out that the samples collected for each Ti are i.i.d. 
and thus the confidence interval can be rigorously constructed based on the evaluation of P for the 
samples. 



4 Sample reuse factor 

Let rij be the number of simulations required at r^. Define sample reuse factor as follows: 

Nl 



where £{X) denotes the expectation of random variable X. Obviously, T TeU se measures the im- 
provement of efficiency upon the conventional method. We demonstrate that the improvement can 
be significant in most applications. 



Theorem 1 The sample reuse factor T reuse = Ijl — Yl\=2 \r^J > w ^ iere d = n for l p ball B p 
and homogeneous star-shaped bounding set Bjj(r); and 

s c 

<* = 2>(?i) + 2>(Ai) 

i=l j=l 

for spectral norm ball B a (r) with k(.) defined as 

k(X) : = 



2mn if X is a variable in C nxm 
mn if X is a variable in R nxm . 
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See the appendix for proof. For illustration purposes, we choose = b — ( b ^ 1 ) f or { — 1 ? 2, . . . , Z. 

By Theorem [0 T reuse = III — Y\\-o 1 t^t 1 • Figures [[] and [2] show that the improvement 

over the conventional approach is significant when d is not large. These figures also reveal that the 
sample reuse factor does not scale well with the uncertainty dimension. For example, when d > 160, 
the efficiency gained from sample reuse techniques may not be attractive. 




O 20 40 GO 80 100 120 140 1 GO 

d 



Figure 1: Performance improvement (A : I = 200, b = 2a; B : I = 100, b = 2a; C : I = 100, a = 
0; D : I = 20, b = 2a). 
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d 



Figure 2: Performance improvement (A : I = 200, b = 2a; B : I = 100, b = 2a; C : I = 100, a = 
0; D : I = 20, b = 2a). 

5 Illustrative examples 

In this section we demonstrate through examples the power of the sample reuse algorithm in solving 
a wide variety of complicated robustness analysis problems which are intractable in the classical 
deterministic framework. 
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First, we consider an example which has been studied in [TT] by a deterministic approach. The 
system is as shown in Figure [3) 



r + e 




C(s) 




P(s) 


c 


— ► 







Figure 3: Uncertain system. 

The compensator is C(s) = and the plant is P(s) = s ( s+ 4+^2S^&+6+o 3<5 3 ) w ^h parametric 
uncertainty A = ^2, ^3] T . The nominal system is stable. The closed-loop roots of the nominal 
system are 

Z! = -15.9178, z 2 = -1.8309, z 3 = -1.1256 + 7.3234?, z 4 = -1.1256 - 7.32341 

The -Hqo norm of the nominal closed-loop transfer function is 1 1 ^Z^ 1 1 cso = 2.78. The peak value, rise 

time, and settling time of step response of the nominal system are, respectively, Pp eak = 1.47, t° = 
0.185, and t® = 3.175. In all of the following examples, we take I = 100. To guarantee that the 
absolute error of the estimate for the proportion is less than 0.01 with confidence level 99%, we 
choose N = 26,492 based on the well-known Chernoff bound (see [12\ [T9] for "sharper" bounds). 
Since the Chernoff bound is conservative, we also performed a post-experimental evaluation of the 
estimates by constructing confidence intervals with confidence level 99% based on Clopper-Pearson's 
method [TP] , 

Figured] is the robustness degradation curve for robust stability over uncertainty set Boo(r) := 
{A : 1 1 Alloc < r }- It demonstrates that a significant enhancement of the robustness margin can be 
achieved at the price of a small risk. 

Figure [5] is the robustness degradation curve, with the robustness requirement P defined as 
stability and norm < 170% HT ^, and the uncertainty set defined as the ellipsoid ^(r) := 
{A : ||A|| 2 < r}. 

Figure [6J is the robustness degradation curve with the robustness requirement P defined as In- 
stability with the domain of poles defined as: real part < —1.5, or it falls within one of the two disks 
centered at z% and Z4 with radius 0.3. The uncertainty set is defined as the polytope 

B H {r) := |rA+(l-r) ^ : A G convjA 1 , A 2 , A 3 , A 4 } j , 

where "conv" denotes the convex hull of A* = [| sin( 2 ^ ll -7r), | cos(^ 1 -7r), — -^] T for i = 1, 2, 3 and 
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Figure 4: Robustness degradation curve (reuse factor = 41). 




Figure 5: Robustness degradation curve (reuse factor = 43). 
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A 4 = [0, 0, 1] T 



Upper Confidence Li 



-> Variance Vari 




Uncertainty Radii_ 



Figure 6: Robustness degradation curve (reuse factor = 49). 

Figure [7] is the robustness degradation curve for the case where the uncertainty set is £>oo( r ) := 
{A : ||A||oo < r}, the robustness requirement P is: stability, and rise time t r < 135% t° = 0.25, 
settling time t s < 110% t° s = 3.5, and overshoot P peak < 116% P° eak = 1.7. 




0.17 0.18 0.19 0.2 



Figure 7: Robustness degradation curve (reuse factor = 38). 
Finally, we consider the same example in [8] where the class of uncertainty is defined as 

A := {blockdiag[<7i/ 5 , q 2 h, Ai]}, 

where Ai 6 C 4x4 , and I5 denotes the identity matrix of 5 x 5. By Theorem [U we have d = 34. 

Figure [8] shows the robustness degradation curve. An improvement (of efficiency) about fivefold is 
achieved by our algorithm. 
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Figure 8: Robustness degradation curve (reuse factor = 5). 

6 A probabilistic perspective 

In sections 2 and 3, we have derived our sample reuse algorithm from the worst-case deterministic 
framework. In this section, we show that the proposed algorithm is also applicable from the per- 
spective of the random nature of uncertainty. In situations where we need to take into account the 
random nature of uncertainty, the pioneering work of Barmish, Lagoa, Tempo, Bai, and Fu [21 [1] 
allows our sample reuse algorithm to be applied to solve efficiently a wide variety of robustness 
problems. The following theorem plays an important role. 

Theorem 2 (see [2]) Suppose that uncertainty A is a random variable with a density function /(A) 
which depends only on £(A) and is nonincreasing with respect to £(A). Then 

Pr{The associated system G(A) guarantees P | A € £>(r)} > inf P(p). 

0<p<r 

Remark 1 . A remarkable fact of Theorem [2] is that no assumption needs to be imposed on the 
robustness requirement P. The assumption in Theorem [2] is roughly interpreted to mean that the 
probability measure of the uncertainty is radially symmetrical with respect to the nominal value. In 
many applications, small perturbations are more likely than large perturbations, and the uncertainty 
is sufficiently unstructured so as to be treated equally likely in the surface of B(r) [2]. 

Remark 2. It should be noted that Theorem [2] applies to a homogeneous star-shaped bounding 
set, l p ball, and spectral norm ball. We introduce in Theorem a conditional probability based on 
the following reason: It does not seem logical to treat the uncertainty as different bounded random 
variables. For example, if the uncertainty possesses a certain distribution over B(ri), it would be a 
contradiction that the uncertainty possesses another distribution over B{r2) for r2 > r\. In fact, if 
the uncertainty is of random nature, then the associated distribution is unique. 

Based on Theorem [2J we can apply the sample reuse algorithm to estimate P(r) for r E [0,6], 
from which we can construct the lower bounds for Pr{G(A) guarantees 
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P I AG B(r)}. 



7 Conclusion 

We develop a fast algorithm for computing the robustness degradation function which overcomes the 
computational complexity and conservatism issue of the deterministic worst-case methods. We also 
demonstrate that our algorithm can provide efficient solutions for a wide variety of robustness analysis 
problems which are intractable by the deterministic worst-case methods. We derive our algorithm 
from the worst-case deterministic framework and also show that the algorithm is applicable from a 
probabilistic perspective. 



A Proof of Theorem 1 

The following lemma follows essentially from the definition of volume function vol(.). 

Lemma 1 Let X = {A £ C nxm : a(A) < r} and Y = {A £ C rnxn : a (A) < r}. Let Z = {A £ 
R" xm : a(A) < r} and W = {A £ R mxn : a(A) < r}. Then vol(X) = vol(Y) and vol(Z) = vol(W). 

Lemma 2 Let m > n. Define spectral norm ball B^'(r) = {A £ C nxm : ct(A) < r} and spectral 
norm ball B^{r) = {A £ R™ xm : a (A) < r}. Then vol(B^(r)) = vol(B^(l))r d with d = 2mn and 
vol(B^(r)) = vol(S^(l))r d with d = mn. 

Proof. By Theorem 1 of [8], we have 

T C ( TT 2(m-n)+l v 1 r / 2 2\2 



f n^ 2(m -" )+1 x [] (af-a 2 k fda 1 da 2 ---da n = l 

Jr>a 1 >a 2 >->a n >0 i=l l<i<k<n 



vol(^(r)) 

with To = rw — r- J rm rrr- Performing a change of variables with X{ = — for i = 1, . . . ,n, we have 



2mn"£ 



vol(S£(r)) 
Thus 



"J x 2 ( m x ^2 _ x 2-j2 dxidX2 • • • dx n = 1. 

->x 1 >x 2 >->x n >0 i=l Ki<k<n 



/n 
JJ x 2 (m -n) + l x Yl (x*-xl) 2 dx 1 dx 2 .--dx n , 
->X!>X2>->x n >0 i=1 l<i<k<n 

and 

vol(s£(r)) = vol(S(l)) r 2mn . 
Similarly, by Theorem 2 of [8], we can show that vol(B^(r)) = vol(i3^(l))r d with d = mn. □ 
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Lemma 3 vol(£>(r)) = vol(B(l))r d where d = n for L ball B p {r) and homogeneous star-shaped 
bounding set Bh{t)', and d = Y2t=i K (li) + Y2j=i K (A?) f or spectral norm ball B a (r). 

Proof. The truth is obvious for cases of an l p ball and homogeneous star-shaped bounding set. To 
prove the lemma for the case of a spectral norm ball, we need to apply Lemmas Q] and [2j □ 

Lemma 4 For i = 1, . . . , I — 1, 



Proof. Let q ,q , . . . , q nj be the samples generated on ry. For I = 1, . . . , nu, define binomial random 



variable Xji+i such that 



1 if q l fall in B(r i+1 ) 
otherwise. 



By the rule of the sample reuse algorithm, 



i aj 

j=i i=i 



Thus for i = 1, — 1, 

i 

£[n i+1 ]=N-J2 £ 



£4 



i+l 



Z=l 



N 



E E E^ i n j = n ] Pr K = n ^ 



j=l n£f2 n . Z=l 



where Q nj denotes the sample space of rij. Since q l is a random variable with uniform distribution 
over B(ri), it follows from Lemma [3] that 



S[X l j>i+1 | n, = n] 



vol(g(r m )) _ ( r i+ i 
vol(B{ rj )) 



Therefore, 



S[n i+1 ] = E n 

7=1 neiin, 



r i+ i 



Pr{n, = re} 



J 7 ra6f2„ 



re Pr{rij = n} 



"-S(t)'' w - 



□ 
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Lemma 5 For i = 2, • ■ ■ , I, 



£ [m] =N-N 



n-i 



Proof. We use induction. Obviously, 



£ [ ni ] = N. 



By Lemma HI we get 



Suppose it is true that 



£ [n 2 ] = N — N 



Then 



i=i 

i 

£ 



N-Nl 



N 



N 



r.i-1 



It follows from Lemma H] that 



3=1 v 3 



£[n i+1 ]=N-J2 l^) £[n j )=N-N 



The proof of Lemma [5] is thus completed by induction. 



□ 



Now we are in the position to prove Theorem [TJ By Lemmas U] and [5l we have 



£ 



,i=l 



i=2 



N — N 



Nl- N 



i=2 



r; 



Therefore, 



Nl 



I 



£[Ei=i n *] i-Y l Jj±- 

and thus the proof of Theorem [T] is completed. 
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